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Abstract 

Ph ' We recall the PFS model constructed for the modeling of unsteady mixed flows in closed water pipes 

■^r , where transition points between the free surface and pressurized flow are treated as a free boundary 

associated to a discontinuity of the gradient of pressure. Then we present a numerical kinetic scheme for 
the computations of unsteady mixed flows in closed water pipes. This kinetic method that we call FKA 
C^ ■ for "Full Kinetic Approach" is an easy and mathematicaUy elegant way to deal with multiple transition 

points when the changes of state between free surface and pressurized flow occur. We use two approaches 
namely the "ghost waves approach" and the "Full Kinetic Approach" to treat these transition points. 
We show that this kinetic numerical scheme has the following properties: it is wet area conservative, 
under a CFL condition it preserves the wet area positive, it treats "naturally" the drying and flooding 
area and most of all it preserves every stationary flow. Finally numerical experiment versus laboratory 
experiment is presented and the scheme produces results that are in a very good agreement. We also 
present a numerical experiment when flooding and drying flows may occur and flnally make a numerical 
study of the order of the kinetic method. 
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Notations concerning geometrical quantities 

k^ , 9{x) angle of the inclination of the main pipe axis z — Z{x) at position x 

^ ' Z(t, x) dynamic slope 

5t 1 il,{x) cross-section area of the pipe orthogonal to the axis z — Z{x) 

S{x) area of il{x) 

R{x) radius of the cross-section il{x) 

(j{x, z) width of the cross-section 51(a;) at altitude z 
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Notations concerning the PFS model 

p(t, X, y, z) pressure 

Pq density of the water at atmospheric pressure po 

p{t,x,y, z) density of the water at the current pressure 

'pit, x) p(t, x) = , . / p(t, X, y, z) dy dz is the mean value of p over VL{x) (press, flows) 

c sonic speed 

S{t,x) "physical" wet area i.e. part of the cross-section area in contact with water (equal to 

S{x) if the flow is pressurized) 

p(t,x) 
A{t,x) A{t,x) — '■ — S{t,x) is the equivalent wet area 

u{t, x) velocity 

Q{t,x) Q{t,x) = A{t,x)u{t,x) is the discharge 

E state indicator. E = ii the flow is free surface, E — 1 otherwise 

'H{S) the .^-coordinate of the water level equal to HiS) — h{t, x) if the state is free surface, 

R{x) otherwise 

p{x, A, E) mean pressure over J7 

Ks > Strickler coefficient depending on the material 

Pm (A) wet perimeter of A (length of the part of the channel section in contact with the water) 

RhiA) Rh{A) = is hydraulic radius 

Pm[A) 

Bold characters are used for vectors, except for S and for Z, the dynamic slope defined later. 

1 Introduction 

The presented work takes place in a more general framework: the modeling of unsteady mixed flows in any 
kind of closed water pipes taking into account the cavitation problem and air entrapment. We are interested 
in flows occurring in closed pipes with non uniform sections, where some parts of the flow can be free surface 
(it means that only a part of the pipe is filled) and other parts are pressurized (it means that the pipe is 
full) . The transition phenomenon between the two types of flows occurs in many situations such as storm 
sewers, waste or supply pipes in hydroelectric installations. It can be induced by sudden changes in the 
boundary conditions as failure pumping. During this process, the pressure can reach severe values and may 
cause damages. The simulation of such a phenomenon is thus a major challenge and a great amount of works 
was devoted to it these last years (see [9, 11, 22, 25], and references therein). 

The classical shallow water equations are commonly used to describe free surface flows in open channels. 
They are also used in the study of mixed flows using the Preissman slot artefact (see for example [9, 25]). 
However, this technic does not take into account the subatmospheric pressurized flows (viewed as a free 
surface flow) which occur during a water hammer. In recent works, [16, 17, 18], a model for mixed flows in 
closed water pipes has been developed at University of Liege, where they use the artifact of the Preismann 
slot for supatmospheric pressurized flow and by introducing the concept of "negative Preissman slot" for 
subatmospheric pressurized flow. They proposed also a numerical scheme to compute the stationary flow, 
as well as the unsteady flow. 

On the other hand the Allievi equations, commonly used to describe pressurized flows, are written in a 
non-conservative form which is not well adapted to a natural coupling with the shallow water equations. 
A model for the unsteady mixed water flows in closed pipes, the PFS model, and a flnite volume discretisation 
have been proposed by the authors in [2] and its mathematical derivation from the Euler incompressible 
equations (for the free surface part of the flow) and from the Euler isentropic compressible equations (for 
the pressurized part of the flow) is proposed in [5]. This model and the flnite volume scheme extend the 
model studied by two of the authors for uniform pipes [6] . In [8] two of the authors has constructed a kinetic 
numerical scheme to compute pressurized flows in uniform pipes. For the case of a non uniform closed pipe 
and for pressurized flow, the authors has extended the previous kinetic numerical scheme with reflections. 



see [1] . Let us also mention that the construction of a kinetic numerical scheme with a correct treatment of 
all the source terms has been published recently [4]. 

The paper is organized as follows. In the second section, we recall the PFS model and focus on the continuous 
flux whose gradient is discontinuous at the interface between free surface and pressurized flow. The source 
terms are also highlight: the conservative ones, the non conservative ones and the source term which is neither 
conservative nor conservative. We use the definition of the DLM theory [19] to define the non-conservative 
products. The particular case of the friction term which is neither conservative nor non-conservative will be 
upwinded using the notion of the dynamic slope, already introduced by the authors in recent works [2, 4]. 
We state in this section the theoretical properties of the system that must be preserved by the numerical 
scheme. 

Section 3 is devoted to the kinetic interpretation of the PFS model thanks to the classical kinetic interpre- 
tation of the system (see [21] for instance). 

In section 4, we construct the kinetic scheme for the PFS model. Firstly, we use the same kinetic scheme 
with reflections that we have constructed in [8, 1, 4] to treat the part of the flow where no transition points 
are present. Then we treat the transition points by two ways: 

• as in [2] , the "ghost waves approach" is used. We make an assumption on the speed of the discontinuity 
between free surface and pressurized flow, to compute the macroscopic states at the right hand side and 
the left hand side of this discontinuity. For this sake, we treat the transition points at the macroscopic 
level. 

• a new approach that we called the "Full Kinetic Approach" is then used to treat these transition 
points. We stay at the microscopic level to build the macroscopic states at the right hand side and the 
left hand side of this discontinuity. 

The particular treatment of the boundaries of the pipes is treated and finally, we construct a well-balanced 
correction to the kinetic numerical scheme that will preserve every stationary fiow. 

In the last section, we present numerical experiments: the first one is the so-called Wiggert's test where we 
have experimental data to compare with. A very good agreement is shown. Then we perform a code to 
code comparison for pressurized fiow: we compare the results of the belier code used by the engineers of 
Electricite de France, Centre d'Ingenierie Hydraulique, Chambery, to compute a numerical solution of the 
Allievi equations by the characteristics method with the one we implemented, called FlowMix, for the same 
engineers for the computation of mixed fiows. We will also test the robustness of the code FlowMix on a 
drying and hooding fiow. The finite volume version of the method we have presented in [2] could not treat 
this type of fiow unless by the introduction of a cut-off function that will produce a lack of conservation of 
mass. Finally we perform a numerical study of the order of the method on a unsteady mixed flow (computed 
by the VFRoe solver that we have constructed and validated in [2]) which will converge to a steady mixed 
flow. 

For the sake of simplicity, we do not deal with the deformation of the domain induced by the change of 
pressure. We will consider only an infinitely rigid pipe (see [7] for unsteady pressurized fiows in deformable 
closed pipe). 

2 A model for unsteady water flows in closed water pipe 

The PFS model (see [2, 5, 12]) is a mixed model of a pressurized (compressible) and free surface (incom- 
pressible) fiow in a one dimensional rigid pipe with variable cross-section. The pressurized parts of the fiow 
correspond to a full pipe whereas the section is not completely filled for the free surface fiow. The free 
surface part of the model is derived by writing the 3D Euler incompressible equations and by averaging 
over orthogonal sections to the privileged axis of the fiow. In the same spirit, by writing the Euler isen- 

tropic and compressible equations with the linearized pressure law p(t, x, y, z) = po H — ^ifii^i ^i V^ ^) ~ Po), 

p(<, x) 
we obtain a Saint- Venant like system of equations in the "FS-equivalent" variable A{t,x) ~ — ^ — Six), 

Pa 



Q(t,x) = A{t,x)u{t,x) which takes mto account the compressible effects (for a detailed derivation, see 
[2, 5, 12]). These variables are suitable to study mixed flows by setting: 



A{t,x) 



Po 



S{t,x), Q{t,x) — A{t,x)u{t,x), 



where S is the physical wet area, i.e. the part of the cross-section area in contact with water. 
In order to deal with the transition points (that is, when a change of state occurs), we introduce a state 
indicator variable E which is equal to 1 if the state is pressurized and to if the state is free surface. 
Notice that S is {A, E) dependent via the relations: 



S = S{A,E) = 



S if E=l, 
A if E = 0. 



The pressure law is given by a mixed "hydrostatic" (for the free surface part of the flow) and "acoustic" 
type (for the pressurized part of the flow) as follows: 



p{x,A,E) = c^ {A-S)+ gh {x,S)cos0 



(1) 



where g is the gravity constant, c the sonic speed of the water (assumed to be constant) and 9 the inclination 
of the pipe. The term Ji is the classical hydrostatic pressure: 

f«(S) 



Iiix,S) = / {H{S)-z)adz 

J-R 



where a{x,z) is the width of the cross-section, R — R{x) the radius of the cross-section and HiS) is the 
z-coordinate of the free surface over the main axis Z{x) (see figure 1 and figure 2). 



= R{x) 



C:z = Z(x) 




Figure 1: Geometric characteristics of the domain: 
free surface and pressurized flow. 




Figure 2: Cross-section il. 



The pressure defined by Equation (1) is continuous throughout the transition points and we define the PFS 
model by: 

f dt{A) + d4Q) = 

(2) 



dtiQ) + d^{=^+pix,A,E)] = -gAZ' + Pr{x,A,E) 



-G{x,A,E) 
~K{x,A,E) 



Q\Qi 
A 



where z = Z{x) is the altitude of the main pipe axis. The terms Pr, G and K denote respectively the 
pressure source term, a curvature term and the friction: 



Pr{x,A,E) = c^ 



A 



1 5' + 9/2(2;, 5') cos 6*, 



G{x,A,E) = gAZ{x,S)=gA{H{S)-h{x,S)/S){cos9Y, 

where we have used the notation /' to denote the derivative with respect to the space variable x of any 
function /(x). The term I2 is the hydrostatic pressure source term defined by: 



l2{x,S)= / {n{S) - z)d^a dz . 

J-R 



The term Kg > is the Strickler coefficient depending on the material and Rh{S) is the hydraulic radius. 
Remark 2.1. The unknown state vector is noted U = (A, Q) and the flux vector F by: 



Fix,U,E)^ [A,^+pix,A,E) 



We simply note (when no ambiguity is possible): 



Q' 



F2{A,Q) = ^+p{x,A,E), 



(3) 



the second component of the preceding flux. 

As it was pointed out in [2, Remark 4-^]^ t^^ flux is continuous through the change of state of the flow 

whereas its derivative with respect to A is discontinuous, due to the jump of the sound speed, see figure 3. 
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Figure 3: Pressure law and sound speed in the case of a rectangular pipe. 
Trajectory (1) corresponds to a pressurization. Trajectory (2) depends on the state of the flow around 



Identification of tiie source terms. 

fn order to write the kinetic formulation of the PFS equations, we have to factorize by A the right hand 
side of System (2). Then, the source terms reads as follows: 



gZ' is a conservative term. 



..^A^S\ ic-[^)S' if E = X 

AS ) \0 if^ = 

/2(x, S)cos9 , 



is a non-conservative product. 



9- 



A 



is a non-conservative product (since I2 could be written as 7(0;, S)S' for some function 
7 specific to geometry of the pipe) . 

• g ('H(S') — Ii{x, S)/S) cos 6*' is a non-conservative product. 

• K{x, A, E) — -pr- is neither conservative nor non-conservative. 

Moreover, all the term said to be non-conservative product are genuinely non-conservative product since 
they do not admit an exact differential form. 

Remark 2.2. We introduce the term Z = [ Z + / K{x, A)u\u\ dx j which is called the dynamic slope 

since it is time and space variable dependent contrary to the static slope Z which is only x-dependent. 
System (2) has the following properties: 
Theorem 2.1. 

1. System (2) is strictly hyperbolic on {A(t,x) > 0} . 

2. For smooth solutions, the mean velocity u = Q / A satisfies: 



9tu + 9x ( Y + c^ HA/S) + gH{S) cos 6* + ,gZ ) = 0. 



(4) 



The quantity ^{A,Q,cos9, Z,E) — ~^ + c \n{A/S) + gH{S) cos 6 + g Z is called the total head. 



3. The still water steady state, for u — Q, reads: 

c^ln{A/S) +gniS) cos 9 + gZ ^cte (5) 

for some constant cte. 
4- System (2) admits a mathematical entropy: 

n2 _ 

£{A, Q,E) = -^ + c^A HA/S) + c^S* + gAZ{x, S) cos 9 + gAZ 

which satisfies the entropy relation for smooth solutions 

dtS + d., {{£ + pix, A, E))u) = . (6) 

Proof. The proof of these assumptions rehes only on algebraic combinations of the two equations forming 
System (2) and is left to the reader. D 

In what follows, when no confusion is possible, the term K{x,A^E) will be noted simply K{x,A) for free 
surface states and K{x, S) for pressurized states. 

Remark 2.3. Equation (5) is the still water steady state equation associated to the PFS equation. Indeed, 
for a pressurized flow (i.e. S — S), when u = Q and A = A{x), the following equations holds: 

c^ ln(^/S') + gR cos + gZ = cte . 

Thus, when S ^ A, Equation (5) provides g'H{A) cos9 + gZ — cte: this equation represents the horizontal 
line for a free surface still water steady states. Moreover, when mixed still water steady states occur, i.e. 
when one part of the flow is pressurized and the other part of the flow is free surface. Equation (5) holds 
again. 

3 The kinetic interpretation of the PFS model 

Recently in [2], we have investigated a class of approximated Godunov scheme for the present PFS model 
in which we show how to obtain by a suitable definition of the convection matrix an exactly well-balanced 
scheme for the still water steady state. We also point out that the upwinding of the source terms into 
the numerical fluxes introduces a stationary wave with a vanishing denominator. We also discuss on the 
possibility to introduce a cut-off function to avoid the division by zero. But, the truncation of the wet area 
A induces a loss of mass which implies a loss of the conservativity property. Moreover, the numerical scheme 
loses accuracy. Therefore, stationary hydraulic jump, drying and flooding area are not accurately computed 
with this kind of numerical scheme. As pointed out in [21], the numerical kinetic scheme are proved to 
satisfy the following stability properties: the water height conservativity, the in cell entropy inequality and 
the conservation of the still water steady state. Unfortunately, it holds only for rectangular geometry. One 
of the main contribution will be to show how to get at least well-balanced scheme for any given geometry, 
the in cell entropy inequality being still an open problem for kinetic scheme in this framework. 
Therefore, the goal of this paper is to construct a Finite Volume- Kinetic (FVK) scheme that will preserve 
every steady states and that will treat "naturally" the drying and flooding area. The big challenge in this 
construction is the fact that the continuous flux has a discontinuous gradient at the interface between free 
surface and pressurized flow, see Remark 2.1. In [2], the finite volume scheme that we have constructed uses 
the "ghost waves approach" to overcome this difficulty. We will see that we can still use this approach for 
the finite volume kinetic scheme but we will prefer to construct a fully kinetic scheme that is a scheme at 
the kinetic level which treats the changes of type of the fiow. 

First of all, let us recall the kinetic formulation of the PFS model based on Perthame's kinetic formulation 
of conservation laws [20] . 

7 



3.1 The mathematical kinetic formulation 

Let X : M — > M be a given real function satisfying the following properties: 



X(w) = x(-w) ^ , / xi^)duj = 1, / w x{^)duj = 1. 

Jm Jr 

It permits to define the density of particles, by a so-called Gibbs equilibrium, 



where b{t, x) ~ b(x, A{t, x), E{t, x)) with 



b{t,x) \ b{t,x) 



(7) 



(8) 



b{x,A,E) = < 



h{x,A) 
g : cos ( 



A 



if E = 0, 



g ^^^'^^K osO + c^ if E = l. 
The Gibbs equilibrium M. is related to the PFS model by the classical macro -microscopic kinetic relations: 

A^ [ M{t,x,Od^, (9) 



Q- / iM{t,x,i)di, 



A 



^Ab{x,A,Ef^ f eM{t,x,0 

Jr 



d^. 



(10) 



(11) 



From the relations (9)-(ll), the nonlinear PFS model can be viewed as a single linear equation involving 
the nonlinear quantity M: 

Theorem 3.1 (Kinetic Formulation of the PFS model). (A,Q) is a strong solution of System (2) if and 
only if Ai satisfies the kinetic transport equation: 



dtM + S, ■ d^M - gcf) d^M = lC{t, x, ^ 
for a collision term )C(t,x,^) which satisfies for (t,x) a.e. 

1 



(12) 



e 



IC{t,x,^)d^ = 0. 



The source terms are defined as: 
with 



{x,W) = B{x,W) ■ d^W 



w ^ (z, s, cose*) 



(13) 
(14) 



c2 fA-S\ -/{x,S)cose -. 



and B 



"^^ g \ AS J A 

l,- '^'''f'''' ,Zix,A) 



,Z{x,S)] ifE^l, 
ifE = 



where l2{x, S) reads 7(2:, S)S' for some function 7 (depending on the geometry of the pipe). 

Proof. The proof relies on very obvious computations since A4 verifies the macro- microscopic kinetic relations 
(9), (10), (11), and from the definition of the source term W and B. D 



Remark 3.1. 

• The kinetic formulation presented in Theorem 3.1 is a (non physical) microscopic description of the 
PFS model. 



• For a circular cross-section pipe, we have: 



1 /H(5)7r . fH{S)\ , aiX,H{S))\ 



2 7rS' V 2 ^ ' \Rix) J 2 y ■ 

4 Construction of the kinetic scheme for the PFS model 

In this section, following the works of [21, 4], we will construct a finite volume kinetic scheme that preserves 
every steady state and that will compute "naturally" flooding and drying zones. The main feature of this 
scheme is the treatment of transition points between free surface and pressurized flows. In a first step, we 
will use the "ghost waves approach" that we have constructed in [2] to treat this difficulty: to this end we 
will go back to the macroscopic level to compute the unknown states {A, Q) at the interface between free 
surface and pressurized flow. 

In a second step, we will construct a fully kinetic scheme to treat the interface between the free surface and 
the pressurized flow: the Gibbs equilibrium on the right hand side and the left hand side of the interface 
between free surface and pressurized flow will be computed by kinetic formulas. To upwind all the source 
terms at the microscopic level, we will use the ideas presented in the recent work of the authors [4] . 
The particular treatment of the boundary conditions will be rapidly exposed using the same approaches. 

4.1 The kinetic scheme without transition points 

In this section, we will treat the parts of the flow that are either free surface or pressurized. Under this 
assumption and based on the kinetic formulation (see Theorem 3.1), we construct easily a Finite Volume 
scheme where the conservative quantities are cell-centered and source terms are included into the numerical 
fluxes by a standard kinetic scheme with reflections [21]. 
To this end, let N € N* , and let us consider the following mesh on [0,L] of a pipe of length L. Cells are 

denoted for every i e [0, A^+1], by to^ = {xi-i/2,Xi^if2), with Xi = and hi — Xi_|-i/2 — a;i-i/2 

the space step. The "fictitious" cells tuq and mj^^i denote the boundary cells and the mesh interfaces located 
at xi/2 and x^v^i^]^ /2 are respectively the upstream and the downstream ends of the pipe (see figure 4). 

2^1/2 < "^V a;Af+i/2 
I \ ^r \ 1 



L 

Figure 4: The space discretisation. 

We also consider a time discretization i" defined by /:"+^ = t" + At" with Ai" the time step. 

We denote [/" = {Af, Q"), u" = — ^, Mf the cell-centered approximation of [/ = {A, Q), u and M on the 

cell rui at time t". We denote by C/q — {Aq, Qg) the upstream and t/Ar+i = (^w+i' Qw+i) ^'^^ downstream 
state vectors. 

For i e [0, A^ -I- 1] , £'i is the state indicator of the cell: ii^i = if in the cell i, the flow is a free surface flow, 
Ei = \ ii in the cell i, the flow is apressurized flow. 



The piecewise constant representation of W defined by (14) is given by, W{t,x) — Wi{t)lmi(x) where 

Wi{t) is defined as Wi{t) — — — / W(t,x)dx for instance. 

Denoting by Wi and M^,+i the left hand side and the right hand side values of W at the cell interface 
2^1+1/21 and using the "straight lines" paths (see [19]) 

*(s, W„ W,+i) = sW,+i + (1 - s)W,, s e [0, 1], 

we define the non-conservative product 0(i,a;i+i/2) by writing: 



t,x,+^/^) = [W]{t)- f B{t,^{s,W,{t),W,+,{t)))ds 
Jo 



Ht,^^+l/2) = [W]it)■ B{t,^is,W,{t),W,+i{t)))ds (15) 

where [W] (i) := Wi+i{t) — Wi{t), is the jump of W{t) across the discontinuity localized at x = a;j+i/2- As 
the first component of S is 1, we recover the classical interfacial upwinding for the conservative term Z. 
Neglecting the collision kernel as in [21, 4] and using the fact that (/) = on the cell rrii (since W is constant 
on rrii), the kinetic transport equation (12) simply reads: 

—M + S, ■ T-A^ = for a; e TOi . (16) 

ot ox 

This equation is a linear transport equation whose explicit discretisation may be done directly by the following 

way. 

Denoting for x £ rrii , f{tn, x,£_) = M^iS.) = -^(^"j Q"i the maxwellian state associated to A" and Q", 

a finite volume discretisation of Equation (16) leads to: 

/r+i(o=Mr(o + ^^ (x-+i(o-xti(^)) (17) 

where the fluxes A^ ■ , i have to take into account the discontinuity of the source term cj) at the cell interface 
Xi+1/2- This is the principle of interfacial source upwind. Indeed, noticing that the fluxes can also be written 

the quantity SMT i — MT i — A^i+i holds for the discrete contribution of the source term in the system 
for negative velocities ^ < due to the upwinding of the source term. Thus SA47 i has to vanish for 

1+ 2 

positive velocity ^ > 0, as proposed by the choice of the interface fluxes below. Let us now detail our 
choice for the fluxes A(r^ i at the interface. It can be justified by using a generalized characteristic method 

for Equation (12) (without the collision kernel) but we give instead a presentation based on some physical 
energetic balance. The details of the construction of these fluxes by the general characteristics method (see 
[10, Definition 2.1]) is done in [12, Chapter 2]. 

In order to take into account the neighboring cells by means of a natural interpretation of the microscopic 
features of the system, we formulate a peculiar discretisation for the fluxes in (17), computed by the following 
upwindcd formulas: 

positive transmission reflection 

-^^+1/2(0 = l{?>0}A^r(O +'Mi<0.i^-2gA^l,^^<0}M7{~d 

+ Iu<o,c2-2<,A0" ,,2>o}-^^+i (-\A'-2ffA0r+i/2 



negative transmission 
negative transmission rcflccti( 



(18) 



Mt,/,io = iu<o}xr+i(o +i{c>o,c2+2,Ac,,/2<o}^r+i(-o 



positive transmission 
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The term A0"j_j^ ,2 in (18) is the upwinded source term (13). It also plays the role of the potential barrier: 
the term ^^ ± 2glS.(f)'^,^,^ is the jump condition for a particle with a kinetic speed ^ which is necessary to 

• be reflected: this means that the particle has not enough kinetic energy ^^/2 to overpass the potential 
barrier (reflection term in (18)), 

• overpass the potential barrier with a positive speed (positive transmission term in (18)), 

• overpass the potential barrier with a negative speed (negative transmission term in (18)). 

Taking an approximation of the non-conservative product <f> defined by Equation (15), the potential barrier 
^'^r+1/2 ^^^ ^^^ following expression: 

AC+1/2 = W] (tn) ■ B (^t,,,^ Q, W-.(i„), W.+l(i. 

which is an approximation of Acj) by the midpoint quadrature formula. 

Since we neglected the collision term, it is clear that f'"^-^ computed by the discretised kinetic equation (17) 
is no more a Gibbs equilibrium. Therefore, to recover the macroscopic variables A and Q, according to the 
identities (9)-(10), we set: 

In fact at each time step, we projected /"(f) on 7V{"(f), which is a way to perform all collisions at once and 

to recover a Gibbs equilibrium without computing it. 

Now, we can integrate the discretised kinetic equation (17) against 1 and f to obtain the macroscopic kinetic 

scheme: 

The numerical fluxes are thus defined by the kinetic fluxes as follows: 



e 



i+i 



Computing the macroscopic state U by Equations (19)-(20) is not easy if the function x verifying the 
properties (7) is not compactly supported. 

Moreover, as we shall see in the next proposition, a CFL conditions is needed to obtain a scheme that 
preserves the positivity of the wet area. 

Proposition 4.1. Let x be a compactly supported function verifying (7) and note [—M, M] its support. The 
kinetic scheme (19)-(20) has the following properties: 

1. The kinetic scheme is a wet area conservative scheme, 

2. Assume the following CFL condition 

At"max(K| + Af 6^ ^^ max/ii (21) 

i i 

holds. Then the kinetic scheme keeps the wet area A positive i.e: 

if for every i e [0, N + 1] , A" ^ then, for every ie [0, iV + 1] , A" ^ 0. 

3. The kinetic scheme treats "naturally" flooding and drying zones. 
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Proof. We will adapt the proof of [21] to show the three properties that verify the kinetic scheme. 
1. Let us denote the first component of the discrete fluxes (20) {Fa)j^,i' 

An easy computation, using the change of variable /i = |^p — 2gA(/)" i , allows us to show that: 



2. Suppose that for every i G [0, A^ + 1], Af > 0. Let us note (,± = max(0,±^) and a 
Equation (17), we get the following inequalities: 



Ai" 
max, hj 



From 






+iu^-23A0,_i/,>o}>'r+i [-y^^ - 2.9A4+1/2 

Since the function x is compactly supported if |^ — u"| ^ Mb^ then Mf{£,) — 0. Thus 

/r+'(e)55 0ifie-<i^M6r, 

as a sum of non negative terms. 

On the other hand, for |^ - <| ^ M6", using the CFL condition < ct|<^| < 1, for ah i, /f +^ ^ since it is 

a convex combination of non negative terms. 

Finally we have Vi G [0,iV + 1] , /" > 0. Since Vi e [0, A^ + 1] , A^'+^ = / f"^\Od^^ we finally get 

V^e [0,Ar + l], A^+i>0. 

3. Suppose Af = 0. Of course, in the mesh i the flow is a free surface flow. So using the definition of M, 

and the fact that the function x is compactly supported, the only term that may cause problem is — r-. 

b(t,xj 



But since 



we get 



A 



A 



b{t,x) y gIi{x,A)cos9 



when A ^ 0, 



lim 



A 



A-fO b(t,x] 

A>0 ^ ' ' 



= 0, 



thus A4"{^) — 0. This is the reason why we say that the kinetic scheme treats "naturally" the drying and 
flooding zones. D 

4.2 The kinetic scheme with transition points 

The transition points are characterized by the points of the pipe where the flow is not of the same type 
on the left hand side and the right hand side of these points. More precisely, they are characterized by 
their localization in the pipe and the speed of propagation w of the change of the sate. We will assume 
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that it exists at most a finite set of transition points and we will consider that these transition points are 
located at the interface between two cells of the mesh, i.e. a point a;i_|_i/2 of the mesh is a transition point 
if Ei 7^ Ei^i. The speed of propagation of the interface defines a discontinuity line x = w t and let us 
introduce U^ = {A^ ,Q^) and U^ = (A+,(5+) the (unknown) states respectively on the left hand side 
and on the right hand side of this line (see figure 5). The speed w is related to the unknowns U by the 
Rankine-Hugoniot condition on the mass equation 



and thus w is given by: 



dtA + d.,Q^Q 



Q+-Q- 



A+ -A- 



^.+1/2 ■ 


^M 


+ 

i+l/2 


u- 

/ 


w = 

/ 

U+ 


Free surface 




"""^ 






/ E, = 


/ 


^.r 


\ 

Press. 





Ui 



H/2 U,- 



Q+ -Q- 

A+ - A- 



Figure 5: Rankine-Hugoniot condition through the line x = wt. 
Free surface state propagating downstream. 



According to the left U and right unknowns U at the interface Xij^i/2 and the sign of the speed w, we 
have to deal with four cases: 

• pressurized state propagating downstream, 

• pressurized state propagating upstream, 

• free surface state propagating downstream, 

• free surface state propagating upstream. 

Assume that U are given then the kinetic scheme in the case of transition points reads: 



where the numerical fluxes are computed by Equation (20) and the numerical microscopic interface quantities 

mulas (18) and the sign of speed w a: 



•^i±i/2 ^^^ obtained according to the formulas (18) and the sign of speed w as follows: 



■^m/2 ^ M-^,^^{tM+,Mf+,) if w<0 



(22) 



M+{^,M-,M- 



if w > 



"^+1/2 ]^ Mt+,^^{^,M+,Mf+,) if u;<0 
where A4^ are the Gibbs equilibrium associated to U according to the formula (8). 
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A first way to compute the unknowns U called "the ghost waves approach" is to go back to the macroscopic 
level and to solve a linearized Riemann problem with discontinuous convection matrix. 
The second way called "Full Kinetic Approach" computes the states A4^ at the microscopic level and the 
state U are recovered by the relations (9)- (10). 

Only two of four cases are considered since we have two couples of "twin cases" : pressurized state is prop- 
agating downstream (or upstream) as shown in figure 6 and free surface state propagating downstream (or 
upstream) as shown in figure 7. 

4.2.1 The ghost v^raves approach 

In order to specify the unknowns U , we have to define four equations. To this end, the ghost waves 
approach is a way to obtain a system of four equations related to the PFS system. Adding the equations 
dtZ = 0, dtCOsO = and dtS = 0, the PFS model (without friction since it is taken into account in the 
term Z in the kinetic formulation (12), (13)) can be written under a non-conservative form with the variable 
V = {Z,cose,S,A,QY: 

dtV + D{V) d^V = 

with D the convection matrix defined by 



W{V) 



where ^(V) = gSOsV-iS) cos0 — c^(^)"^ and u = Q/A denotes the speed of the water. c{V) is then c for 

the pressurized flow or < /g cos 9 for the free surface flow. 

Remark 4.1. Let us remark that, since dxIi{x,A) = l2{x,A) + dAlii^, A)dxA, the pressure source term I2 
does not appear. 

We assume that the propagation of the interface (pressurized-free surface or free surface-pressurized) has 
a constant speed w during a time step. The ghost waves approach consists to solve a linearized Riemann 
problem in each zone (see [2]): 



/ 











\ 












































1 


^ gA 


gAHiS) 


*(y) 


c'{V)-u^ 


2u J 



dtV + Dd.xV = 









V 




if 
if 


X < wt 
X > wt 



where D — D{V). The term V is some approximation depending on the left Vi and the right Vr state. 
The half line x — wt is then the discontinuity line of D. Both states Ut and U (resp. Ui+i and U^) 
correspond to the same type of flow. Thus it makes sense to define the averaged matrices in each zone as 
follows: 

• for X < wt, we set Di = D[Vi) for some approximation Vi which connects the state Vi and V . 

• for X > wt, we set Di+i = D(Vi+i) for some approximation Vi+i which connects the state V and 

Then we formally solve two Riemann problems and use the Rankine-Hugoniot jump conditions through the 
line X — wt which writes: 

Q+-0- = w{A+-A-) (23) 

F2{A+,Q+)^F2{A-,Q-) = w{Q+-Q-) (24) 
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with F2 defined by (3). In what follows, all quantities of the form Vi (resp. Wi+i) define some approximation 
which connects the state Vi and v~ (resp. u"*" and Wi+i). States can be connected, for instance, by the mean 
value of each state. 

Pressurized state propagating downstream: it is the case when on the left hand side of the line 
^ = wt, we have a pressurized flow while on the right hand side we have a free surface flow and the speed 
w of the transition point is positive. Following Song [24] (see also [13]), an equivalent stationary hydraulic 
jump must occur from a supercritical to a subcritical condition and thus the characteristics speed satisfies 
the inequalities: ^ 

Ui+l + c{V)i+i < w <u,+c. 



£. = Ui — c 




^ = Ul 



.. 5 = Ui + c 



Figure 6: Pressurized state propagating downstream. 

Therefore, only the characteristic lines drawn with solid lines are taken into account. Indeed they are related 
to incoming waves with respect to the corresponding space-time area — c» < ^ < w. Conversely, the dotted 
line ^ = Mj+i — c{V)iJ^i, for instance, related to the free surface zone but drawn in the area of pressurized 
flow is a "ghost wave" and is not considered. Thus U~^ = Ui+i and Ui, U~ are connected through the 
jumps across the characteristics ^ = and S, =Ui — c. Eliminating w in the Rankine-Hugoniot jump relations 
(23)-(24), we get U^ as the solution to the nonlinear system: 



(^2(A,+i,Q,+i)-F2(A-,g-)) = 



Q -Q.,-{A -Ai){u,-c) + 



{Q^+l-Q-? 

{A,+i~A-) 



C + Ui 



= 



(25) 



where tpl^ is the upwinded source term 

Z,+i - Zi+H{S)(cos0,+, - cos^,) + ^{V){S,+i - S,) . 

Free surface state propagating downstream: it is the case when on the left hand side of the line 
S, = wt, we have a free surface flow while on the right hand side we have a pressurized flow and the speed w 
of the transition point is positive. Following again Song [24], the characteristic speed satisfies the inequalities: 

Ui + c{V)i < W < Ui+i + c 
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5 = 



^ = w 



? = «»+i - c 




? = '«<+! + t 



Figure 7: Free surface state propagating downstream. 

There are two incoming characteristic Unes with respect to the free surface area — cxd < ^ < u; (actually 
three with ^ = 0) and they can connect the given left state [/,; with any arbitrary free surface state. Thus 
only one characteristic line (^ = Ui+i + c) gives any information (it is Equation (27) above) as an incoming 
characteristic line with respect to the pressurized zone w < S^ < +cx). From the jump relations through 
the characteristic ^ = 0, and after the elimination of w in the Rankine-Hugoniot jump relations (23), (24) 
we get another equation, namely Equation (28) above. It remains to close the system of four unknowns 
{A~ , Q^ , A^ , Q^)- Thus, from Equation (4), we use the jump relation across the transition point (with 



speed w) for the total head ^{A, Q, cos 9, Z, E) = 



'In 



g TilA) cos 6 + g Z, which writes: 



$H 



$ — w (u^ — u ) 



(26) 



Finally, we use the relation: 



W = Wpred with Wpred 



Q 



i+l 



A,+i - A, ' 
which is the predicted speed of the discontinuity. We have then to solve the nonlinear system: 

{Q^+l - Q+) - [A^+i - A+) {u,+i + c) 



(27) 



{Q+ - Q-) iQ^+l - Q^) = {A+i - A,) {F2{A+ ,Q+) - F2{A- ,Q-)) 

^^^^' ' c''\n{A+) + gcos0HiA+) ^ }f J^ ^ c^ln{A-) - gcos9H{A-) 

2 [A ) 



2{A+y 



Q, 



+1 



Q^ fQ+ Q- 



A.,+i-A, \A+ A 
{Q^+l - Q,:) iA+ - A-) = (Q+ - Q-) (A,+i - A,) 



(28) 



(29) 
(30) 



4.2.2 The Full Kinetic Approach 

As we will show in the numerical experiments part of this article, the "ghost waves approach" method seems 
to produce very good results in agreement with experimental data. But, as we have shown in the previous 
section, to solve the interface problem between the free surface and the pressurized part of the flow, we have 
expressed mathematical relations between the macroscopic unknowns A and Q whereas in the free surface 
part as well as the pressurized part of the flow, we only deal with the microscopic quantities A^" and we 
have constructed kinetic fluxes at the interface that take into account every source terms. 
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At this point, we proposed to solve the interface problem at the microscopic level. This method is based on 

the generalized characteristics method applied to the transport equation (12), see [10, 12]. 

As in the previous method we will only consider two generic cases: pressurized state propagating downstream, 

as shown in figure 8(a), and free surface state propagating downstream, as shown in figure 8(b). For these 

two cases, we have to determine the Gibbs equilibrium M corresponding to the macroscopic states U on 

both sides of the curve T representing the trajectory of the transition point whose unknown velocity is w. 

On each time step, this trajectory is supposed to be a line and w must be defined such that the Rankine 

Hugoniot relations are satisfied. 

We will now use the Gibbs equilibrium Ai^ instead of A^i+i in the computations of the microscopic fiuxes, 

see formula (22), for the case w > 

The transport equation (12) is used only for the transmission cases since we are only interested in trajectories 

which intersect F. 

At the microscopic level, we obtain: 

V^>u;, M-{0 l{«^+2sA0^^,/,>o} = A^r {^e + 25 A(/.^Yi/2 ) 1u^+2sA0^^,/,>o}, 
and 

y^<w, x+(o = Mr+i(e)- 

These relations give, at the macroscopic level: 



^^^M-{Otie+2,A4>-,,,,>o}d^ = J^°° (^l^M7{^e + 2^^€^2) he+2,A^^^,,,>o}dC 



and 



£(J)A.^(e).e^£(J)^r..(o^e 



(31) 



(32) 



Since we have seen in Proposition 4.1 that in order to preserve the positivity of the wet area A^ the scheme 
needs to verify a CFL condition (thus a compactly supported function x), for the direct computations of 
all integral terms, we have chosen in the industrial code FlowMix (used by the engineers of Electricitc de 
France, Centre d'Ingenierie Hydraulique, Chambery): 



1 



X = 7wsl[-x/3,V3] • (33) 



Let us mention that some preceding relations are obvious since the support of A^(0 is [u ± c{A)y3]. For 

instance for a free surface fiow propagating downstream (or for the twin case of a pressurized flow propagating 

upstream), in the free surface zone, we have to compare in Equations (31) the velocity of the interface w 

with u + c{A)\/3 while in the pressurized zone u + c(A) v3 is very large . The only corrective term is due to 

the slope. 

From now on, we omit the indexes n to simplify the notations and every computations are done with the 

density function x defined by (33). 

Let us define the "effective" boundary in the integral used in Equations (31): 

w' — max I w, 4/25 max(0, — A(/)j^i/2) ) (34) 



w 



(^max(0,w;2 +2(7A0^^^/2), ^2 5 max(0, A0,+i/2)) . (35) 



Equations (31) write: 

A- 



(S--!-) - ^ (v^*^.' - 2g A0,+i/2 - ^-ff - 2g A0,+i/2) 



c{A 



(36) 



{{n^-ii-Y) = -fr-AK-it) (37) 
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where 7 , S , 7^, Si are defined by: 

7^ = max (w', u" — c(A^)v3J , 6^ — max {w',u^ + c(A^)v3 j 



and 

7i = max f ui", u^ — c(Ai)v3 ) , 5i = max 

In the same way, Equations (32) write: 



\w",u^ + c{Ai)V2, 



(/3^ - a^) = , l"^ . (ft+i - a^+i) 



c(A+) 



c(A 



i+l) 



A+ 
liA+ 



{{P^r~{a+r) - ^^{Pl,~al,) 



c{A 



i+ij 



(38) 
(39) 



where a"*", /3+, ai+i, /3i+i are defined by: 



a = nrm 



and 



(w, u+ - c{A+)Vtj , ^+ = min (w, u+ + c{A+)Vtj 

i+c(A,+i)\/3) 



tti+i = min (w, w,j+i - c{Ai+i)V3] , Pi+i = min ( w, Ui+ 
Taking into account Equation (38), Equation (39) becomes: 

a+ + (5+ = Ui+i + /3i+i 
The unknowns are still w, A^ , Q^ , A+, Q+. 



Pressurized state propagating downstream: let us first remark that in this case, since in the mesh i 
the flow is pressurized, we always have w' < Ui + c{Ai)\/3 so that 7^ — (5^ 7^ 0. We solve then the system 
formed by the Rankine-Hugoniot jump conditions (23)-(24), and Equations (36), (37), and (38) which formed 
a full nonlinear system of 5 equations with 5 unknowns. So this procedure privileges the information coming 
from the zone containing the cell interface Xi-f-i/2- 

Free surface state propagating downstream: again, we solve the nonlinear system formed by Equa- 
tions (23)- (24), (36), (37) and (38) and the experience acquired with the industrial code FlowMix has shown 
that the system is a full nonlinear system of 5 equations with 5 unknowns. 

But since in the mesh i the flow is a free surface one, we may have the critical case where w' > Ui + c{Ai)\/3 
with w' deflned by Equation (34). In this case j~ ^S^ =0 and the system formed by the nonlinear equations 
(23)-(24), (36), (37) and (38) is under determined. In this case, as in the ghost waves approach, we replace 

w by Wpred — — j J- and we solve the system formed by the Rankine-Hugoniot jump conditions (23), 

Ai+i — Ai 

(24) and (26), as in the ghost waves approach, completed by Equations (38) and (39). 
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(a) Pressurized state propagating downstream. 



(b) Free surface state propagating downstream. 



Figure 8: Full Kinetic Approach for the case of a transition point. 



4.3 Boundary conditions 

Let us recall that a;i/2 and a;jv+i/2 are respectively the upstream and the downstream ends of the pipe. At 

this stage, we have computed all the "interior" states at time t"+-^, that is {U"~^^)i^i^N are computed, and 

also the state of the "interior" cells {E^~^^)i^i^N by the method presented in section 4.4. 

The upstream state Uq corresponds to the mean value of A and Q on the "fictitious" cell too = (a;_i/2, a;i/2) 

(at the left of the upstream boundary of the pipe) and the downstream state U n+i corresponds to the mean 

value of A and Q on the "fictitious" cell toat+i = (2;Ar+i/2,a;jv+3/2) (at the right of the downstream end of 

the pipe). 

Usually, we have to prescribe one boundary condition related to the state vectors C/g and C/^_|_^ (there is 

generically only one incoming characteristics curve for the upstream and only one outgoing characteristics 

for the downstream). For instance, at the upstream end of the pipe, one of the following boundary conditions 

may be prescribed (we omit the index n + 1 for the sake of simplicity) : 

1. the water level is prescribed. So let Hup{t) be a given function of time. Then we have: 

"' ^Aoity 



Vi > 0, — In 
9 



So{t) 



niSoit))coseo 



Hup {t) 



2. the discharge is prescribed. So let Qup{t) be a given function of time. Then we have: 

Vi > 0, Qo{t) = Qupit) . 

3. the total head may be prescribed. So let $np(^) be a given function of time. Then we have: 



Vt>0, 



2Ao{t) 



=2 In 



So{t) 



gUiSoit)) cos 6*0 + gZo = <^upit) 



(40) 



(41) 



(42) 



At the downstream end, similar boundary conditions may be defined. In order to find a complete state 
for the upstream boundary, Uq , and the downstream boundary, V"^,^, we have to define the missing 
equation b^p and hdown respectively. We will present the method at the upstream boundary of the pipe (it 
is easy to adapt it at the downstream boundary). 
At this stage, we have to consider two cases: 

1. the upstream boundary is not a transition point, that is Eq — Ei, 

2. the upstream boundary is a transition point, that is Eq =/= Ei. 

We will treat the first case at the microscopic level while for the second case, we will describe the "Full 
Kinetic Approach" (the "ghost waves approach" for interior cells may be easily adapted). 
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4.3.1 The boundary is not a transition point 

The jump A(j)^,^ is thus computed on the first half mesh. 

Figure 9 represents the case when the upstream boundary of the pipe is not a transition point between free 
surface and pressurized flow. The Gibbs equilibrium Aig must be determined and we recall that we know 
either Aq (if the upstream water level is prescribed by Equation (40)), Qo (if the upstream discharge is 
prescribed by Equation (41)) or a relation between these quantities (if the upstream total head is prescribed 
by Equation (42)). 




-^e-2gAr„: 



Outgoing 
transition point 



tl/2 



Ml 



Figure 9: The upstream boundary of the pipe is not a transition point. 

To obtain another equation denoted bup{AQ, Qo), we will use the information transmitted from the outgoing 
kinetic characteristic (see figure 9). At the macroscopic level, we obtain: 



?0 



MoiOdC 



^ M-,[-^e-2gAcj^-^^ 



dt 



with ^0 = ^\ 1^9 max(0, A(/)",2)- This also writes as: 



c(Ao) 



{5o - 7o) = 
(^0^-70^) = 



c(Ao) 
where 70, 5q, 71, ^1 are defined by: 

70 = min (^^0 



Ax 

c{A,) 

c(Ai) 



{51 -ll) 



52 + 2gA0i/2 



and 



71 



Wo - c(Ao)V3J , 5o ^ min f ^0, "0 
min Ui,""! - c(Ai)a/3J , (5i = min f^i, ui 



3(Ao)V3 
z{Ax)Vi) 



(43) 
(44) 



and ^1 — — a/Si? max(0, — A0",2)- 

So we will use only one of the two equations (43) or (44) as hup- 

Let us mention that in the case when 71 = (5i, Equations (43) and (44) are useless. This happens when: 

• the flow is an incoming supcritical free surface flow at the upstream boundary condition (this could be 
due to a high slope inducing a great |^i|). In this case, we impose the critical flow that is uq — c{Aq). 

• the flow is an outgoing supcritical free surface flow at the upstream boundary condition. The boundary 
conditions are then useless (two outgoing characteristics), so we impose Aq = Ai, Qq = Qi. 

Remark 4.2. We will use 



20 



1. Equation (43) "momentum of order 0" if Qq is prescribed at the upstream end, 

2. Equation (44) "momentum of order 1 " if Aq is prescribed at the upstream end, 

3. Equation (44) if the total head is prescribed. Since Equation (42) couples the two unknowns, we had 
the choice between the two equations (43) and (44) . The experience acquired with the industrial code 
FlowMix makes us to use Equation (44). 



4.3.2 The boundary is a transition point 

We suppose that Eq ^ Ei. In the first step, we apply the procedure described in section 4.2.2: the left state 
is the downstream state Uq, known at time i", and the right state is C/i known at time t""*"^. We determine 
the state U~ at the left of the transition curve. As the two states C/q and U~ represent the same type of 
flow (free surface or pressurized), we apply the preceding method by just replacing A^" by Ai^ in all the 
formula above (see figure 10). 



Mo .. * 


J 




i '"""- 


"■■■■ 






-^e-29Ar,/y 






y^ Incoming 




z 


/ transition point 
*- 



^1/2 Ml 

Figure 10: The upstream boundary of the pipe is a transition point. 



4.4 Updating the state of the mesh E 

To update the state of the mesh (see figure 11), we use a discrete version of the state indicator E equal to 1 
for a pressurized flow and otherwise. Following [6], after the computation of the wet area ^"^^ we predict 
the state of the cell rtii by the following criterion: 

• if E'^ = then: 

if A^'+i < S., then £;f +i = 0, else E'^"^^ = 1, 



• if ^f = 1: 

if A"+i > S, then £;r 



1, else £;" = £;" 



77' n 



Indeed, if A"^^ > Si it is clear that the mesh rtii becomes pressurized, on the other hand if A"^^ < S'i in a 
mesh previously pressurized, we do not know a priori if the new state is free surface (p — po a-nd the value 
of the wet area is less than Si) or pressurized (in depression, with p < pq and the value of the wet area is 
equal to 5*^: see Remark 4.3 and figure 12). 

So far as we do not take into account complex phenomena such that entrapment of air pockets or cavitation 
and keeping in mind that the CFL condition, defined by (21), ensures that a transition point crosses at most 
one mesh at each time step, we postulate that: 

1. if the mesh rui is free surface at time i", its state at time i"+^ is only determined by the value of A"^^ 
and it cannot become in depression. 

2. if the mesh rui is pressurized at time i" and if A"^^ < S'i, it becomes free surface if and only if at least 
one adjacent mesh was free surface at time t". This is exactly the discrete version of the continuous 

— criterion explained in Remark 4.3 and displayed on figure 12. 
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£j.+i ^ 



yes, if A1+^ < S, 




t — tji 



t = t„+i 




Figure 11: Update of the state E" of the mesh m 



Remark 4.3. During a pressurized state, we may have A > S or A < S . The case A > S corresponds to 
an overpressure while the second one is observed for the depression. As said before, the main difficulty to 
detect depression comes from the fact that the criterion A < S corresponds also to a free surface state. To 
dissociate these two cases, we have to know if p — po, that is also -g = 1 or not. On figure 12, we represent 

... . ... „ ....„„, „„.„.»„ „„.;.. ...„. ... .,. „. I .„.„„„. .„ .„ ... „ ..„ 

a physical situation at different time i^, i = 0, . . . , 3. We draw the behavior of the interface speed w in the 
{x, t) -plane and the graph of the function — at fixed tim,e t^. 
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t = to 



t = U 



t = t2 



t = h 




dx 



w = 



Pres§. 
Supaitm 



flow 



dx 



w = 




Subatm.. 

■""a < s "■■■■ 
s = s 

p(x,A,l) <0 



Figure 12: — criterion to detect depression area. 



4.5 A well-balanced correction 

In this section we make a review of the features of the numerical scheme and proposed an improvement 
of the numerical kinetic scheme that wih preserve every stationary state. First of aU, we recall that given 
any compactly supported function x verifying the properties (7), the numerical kinetic scheme is wet area 
conservative, preserves the positivity of the wet area under the suitable CFL condition (21) and finally it 
treats "naturally" the flooding and the drying zones. 
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But as pointed out by Bourdarias et. al in [1], there is only one choice for the function x that will ensure 
the "well-balanced" properties that is the preservation of the still water steady state {u — 0). Unfortunately, 
Bourdarias et. al in [1, Theroem 3.2] showed that for pressurized flow, this choice is: 



x(^) 



1 / UJ 

exp — TT 



2tt 



which is not compactly supported. For the free surface flow in non uniform closed water pipes, Ersoy 

has shown in his PhD thesis that the function x that will solve the still water steady state is not known 

"analytically" and cannot be computed easily. 

This fact conducts us to use the same type of corrections we have constructed for the VFRoe numerical 

scheme: an exactly well-balanced scheme (see [2, Theorem 6.1]). Thus we have to correct the fluxes and 

constructed cell-interface states in order to preserve every steady states. 

For i € [0, A^ -I- 1], we introduce stationary profiles {Af_^,A~^-^^), see figure 13, stationary states Uf_-^ = 

(A^;^,Qi-i_i), and associated source terms {Zf_i, Z~^^) defined by: 



$(t/,+i,cos6i,+i,Z^^i,^i+i) = <P{Ui+i,cosei+i,Z,+i,E,+i) 
$(t/+_i,cos6ii_i,Z+_i,£;i_i) = $([/i_i,cos6l,_i,Z,_i,£;,_i) 



where the couple [Z^j^^^ Z^_-^) is defined by: 

^r+i = 
zU = 



z, 

Zi+1 


if 
if 


A>i ^ A 


z, 

Zi-i 


if 
if 


AU = A 

Ati + A 



(45) 



(46) 




Figure 13: Steady states profiles. 



In [23] the authors construct a "stationary profiles" finite volume scheme to preserve every steady sate. In 
the same spirit, we consider the following numerical scheme: 



u: 



71+1 



u? 



Ar 



F;f^(u,,Ui^,,Z,,Zr^,;A<l^^,/,) - F+_:-^(uU,U,,ZU,Z,-,A<Pti 



/2 



where the numerical fiuxes are computed by: 



^^+h 



S )-^ii(Orfc. 



(47) 



(48) 



with M-;1 - M-;l ([/„ C/-+1, Z„ Z-+i; A^^^^z^) and X+J^ {uU,U,, Z^, Z,; A<p^_^^^) . Let 
tion that we have explicitly written the dependency of the source terms in the fluxes. 

Theorem 4.1. The numerical scheme (47) preserves every discrete steady states. 
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Proof. According to the first equation of System (2) and to Equation (4), a stationary steady state is 
cliaracterized by a constant discharge and a constant total head along the pipe which writes: 

y X e[0,L],Q{x,E)= cte and $(A, g, cos 6*, Z, E) = cte . 

The goal is to prove that if the states, (t^")o<j<Arj_i £it some time i", verify the following discrete equilibrium 
property: 

Vi e [0,iV + l] , Qi = cte and ^{Ai, Q, , cosOt, Zi,Ei) = cte (49) 

then we have: 

Vm >n,(C/r)o<.<A.+i = (C^r)o<.<A.+i • 

Assume, by induction, that the states (C^")o<i<Ar+i verify the discrete equilibrium property (49). We want 
to show that: 

Vto > n, {UT)o<^<N+l = iU7)o<^<N+l ■ 

By construction the numerical fluxes defined by (48) are consistent. Indeed, we have: 
VU,Z FT^^^^iU, U, Z, Z; 0) = F+_y^{U, U, Z, Z; 0) - F{U) . 

Since V i G [0, A^ + 1] , Q" = cte, and the second component of the flux is monotone with respect to A, we 
deduce from Equations (45) and (46) that Af_-^ = Ai, A~j^-y = Ai, Zf_^ — Zi and Z^^^ = Zi which ensure 
from the consistency of the numerical fluxes that: 

and hence (C/""'"^)o<,;<7v+i = (C/r)o<i<JV+i- □ 

5 Numerical experiments 

This section is devoted to a numerical validation of the model and the kinetic numerical scheme for three 
main cases: single point pressurized flow, the Wiggert's test case, a code to code validation for pressurized 
flow in uniform pipe and an example of drying and flooding flow. At the end of the section, we present a 
numerical study for the order of the discretisation. 

5.1 Single point pressurised flow 

In this section, we present numerical results for the case of a single point pressurized flow, namely the test 
proposed by Wiggert [26]. The numerical results are then compared with the experimental ones: a very 
good agreement between them is shown. According to the experimental data, we were able to propose a 
value (or a range of values) for c which seems physically relevant. On the contrary, in the Preissmann slot 
technique ([26, 14]) the value of c is related to an arbitrary value (the width of the slot) and cannot exceed 
practically lOm/s, otherwise the method becomes unstable. The following test case, is due to Wiggert [26]. 
The experimental device (see figure 14) is an horizontal 10 m long closed pipe with width 0.51 m and height 
H = 0.148 m. The Manning number is 1/Kg — 0.012 sjm^l'^ . The initial conditions are a stationary state 
with the discharge Qo = and the water level h^ — 0.128 m,. 

Then a wave coming from the left side causes the closed channel to pressurise. The upstream condition is 
a given hydrograph (j/2 in figure 16), at the downstream end, a step function is imposed: the water level is 
kept constant to Kq — 0.128 rn until the wave reaches the exit. At this time, the level is suddenly increased 
(see 2/3 in figure 16). For the computations, these boundary conditions have been read on Wiggert's article 
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and rebuilt using piecewise polynomial interpolations (figure 15 below). 

Other parameters are: 

Discretisation points : 80, 

Delta X (m) : 0.125, 

CFL : 0.5, 

Simulation time (s) : 18, 

Sound speed {ms~^) : 40. 



10 m 



wood ^ 
y ^ ^ ^ ^ ^ ^ ^ ^ X ^ 



gate 



concrete 



glass 



14,8 cm 



Figure 14: Experimental device (adapted from Wiggert [26]). 




Figure 15: Wiggert's test : upstream hydrograph (up) and downstream water level (down). 



Ha 5-;a 




^ — 40sec-J 



Figure 16: Wiggert : experimental data. yi : upstream hydrograph, yj, : downstream hydrograph. Iia^ hs, 
he, ho : pressure head at 0.5 tti, 3.5 tti, 5.5 m and 9.5 m from the tunnel entrance (location of recording 
instruments) ([26]). 
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Let us define the piezometric head by: 



piezo=Z + 'H+p with 



P = 



(? (P - Po) 
Po.g 



if the flow is pressurised 



p = h the water height if tlie flow is free surface 



In figure 17(a), we present the piezometric line computed at 3.5 m from the tunnel entrance (solid curve). 
Circles represent experimental data read on curve hs , including maxima and minima points of the oscillating 
parts. We can observe a very good agreement with the experimental data even for the oscillations. We 
point out that we did not find in other papers, by authors carrying out the same simulation, a convenient 
numerical reproduction of these oscillations : they do not treat the dynamical aspect of the pressurized 
flow, in particular when using the Preissmann slot technique ([26, 14]). On the other hand, we found in M. 
Fuamba [13] a similar and interesting approach with a non conservative formulation and another numerical 
method (characteristics). 

The value of the sound speed c was taken equal to 40 m/s^ roughly according to the frequency of the 
oscillations observed during the phase of total submersion of the tunnel. This low value can be explained by 
the structure of the tunnel and by bubble flow (see [15, 25] for instance). 

We observe that the front reaches the control point at 3.6 s, in a good agreement with the experimental 
data (less than 0.15 s late). Let us mention that before it reaches the exit (part AB in figure 17(a)) the 
oscillations of the pressure associated with the moving front reflect between upstream and the front itself 
(since the free surface is at constant pressure) where the channel is flooded. Beyond point B the oscillations 
result from the step in the downstream water level and they propagate in the fully pressurized flow (their 
frequency was estimated using the BC part of the experimental curve). 

Figure 17(b) gives the evolution of the front's speed. We observe the same behaviour as in [26, Figure 7]: the 
front quickly attains a maximum speed, decelerates and then slowly accelerates as it approaches the tunnel 
exit. Moreover the values are consistent with those of Wiggert. 
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Figure 17: Numerical results for Wiggert's test. 



5.2 Numerical validation for a pressurized flow in an uniform pipe 

We present now numerical results of a water hammer test. The pipe of circular cross-section of 2 m^ and 
thickness 20 cm is 2000 m long. The altitude of the upstream end of the pipe is 250 m and the angle is 5°. 
The Young modulus is 23 10^ Pa since the pipe is supposed to be built in concrete: thus the sonic speed is 
equal to c = 1414.2 m/s. The total upstream head is 300 m. 
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Other parameters are: 

Discretisation points : 1000, 

Delta X (m) : 2, 

CFL : 1, 

Simulation time (s) : 100, 

Sound speed (ms'^) : 1414.2. 

To ensure that the model and the kinetic numerical method that we propose describe precisely flows in closed 
uniform water pipes, we present a validation of it by comparing numerical results of the proposed model with 
the ones obtained by solving Allievi equations by the method of characteristics with the so-called belier 
code used by the engineers of Electricite de France, Centre d'Ingenierie Hydraulique, Chambery, [27]. 
A first simulation of the water hammer test is done for a fast cut-off of the downstream discharge for a pipe 
whose Strickler coeeficient is Kg — 90: the initial downstream discharge is 10 m'^ /s and we cut the flow in 
5 s. In figure 18, we present a comparison between the results obtained by our kinetic scheme scheme and 
the ones obtained by the belier code at the middle of the pipe: the behavior of the piezometric line and 
the discharge at the middle of the pipe. One can observe that the results for the proposed model and the 
numerical kinetic scheme are in very good agreement with the solution of Allievi equations. 
A second simulation of the water hammer test is done for the same rapid cut-off of the downstream discharge 
but for a frictionless pipe whose Strickler coefficient is Ks = 215, 6310^. In figure 19, we present a comparison 
between the results obtained by our kinetic scheme scheme and the ones obtained by the belier code at the 
middle of the pipe: the behavior of the piezometric line and the discharge at the middle of the pipe. One 
can observe again that the results for the proposed model and the numerical kinetic scheme are in very good 
agreement with the solution of Allievi equations. 
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Figure 18: Piezometric line (left) and discharge (right) at middle of the pipe. 
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Figure 19: Piezometric line (left) and discharge (right) at middle of the frictionless pipe. 
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5.3 Numerical validation for drying and flooding flow 

We present now numerical results for a flow that will be drying and flooding. The frictionless pipe is 
constituted by a pipe of circular cross-section of diameter 2 m and 50 m long with slope 0.003 and another 
pipe of circular cross-section of diameter 2 m and 100 m long with slope 0.05 . The altitude of the upstream 
end of the pipe is 100 m. The upstream and downstream discharge is kept to 0. 
Other parameters are: 

First pipe discretisation points : 100, 

Delta X (m) : 0.5, 

Second pipe discretisation points : 200, 

Delta X (m) : 0.5, 

CFL : 0.9, 

Simulation time (s) : 500, 

Sound speed {ms~^) : 10. 

The initial state is a flow of constant height (1.8 m) on half the first pipe and a dry zone on the rest of the 
pipe, see figure 20. We present the flow at time T = 6 s, see figure 21, where a drying zone is present, at 
time r = 80 s, see figure 22, when the flow has reached the downstream end and is partially pressurized 
and the flow at the final time T = 500 s, see figure 23 where all the water is in the second pipe. This non 
physical test shows that the kinetic numerical scheme treats "naturally" the drying and flooding zone. The 
water height is exactly equal to 0, in the initial condition and at the final time for the dry zones. 
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Figure 20: Piezometric line (left) and discharge (right) at initial condition. 
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Figure 21: Piezometric line (left) and discharge (right) at time T = 6 s. 
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Piezomelric line (m) at time T = 80 s 
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Figure 22: Piczomctric line (left) and discharge (right) at time T = 80 s. 
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Figure 23: Piezometric line (left) and discharge (right) at final time T = 500 s. 

5.4 "Ghost waves approach" versus "FuU Kinetic Approach" 

We want to compare numerically the two approaches on a violent water hammer "numerical" test for a non 

uniform frictionless closed water pipes. 

To this end, the numerical experiment is performed in the case of an expanding 5 m long closed circular 

water pipe with slope. The upstream diameter is 2 m and the downstream diameter is 3.2 m. The altitude 

of the main pipe axis is set io Z = Im. 

At the upstream boundary condition, the piezometric line (increasing linearly from 1 to to 3.2 to in 5 s) is 

prescribed while the downstream discharge is kept constant equal to Om? /s (see figure 24). The simulation 

starts from a still water free surface steady state where the height of the upstream is 1 to, (see figure 24) and 

the discharge is null. 

Other parameters are 

Discretisation points : 100, 

Delta X (to) : 0.05, 

CFL : 0.8, 

Simulation time (s) : 5, 

Sound speed (tos^^) : 20. 

Let us mention that we have already used this numerical test case to compare the kinetic scheme using the 
"ghost waves approach" with the VFRoe scheme presented in [2], see [4, Figure 3]. 

This numerical test intends to reproduce a "sharp" water hammer experiment inducing large oscillation of 
the piezometric level and the discharge as showed in figures 25 and 26. From a numerical point of view, it is 
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a "hard" numerical test. In order to validate numerically this approach and due to the lack of experimental 
data in the case of variable cross section pipes, we compare the result of the presented numerical scheme 
with those obtained by the upwinded VFRoe scheme [2] . Results are represented in figures 25 and 26 where 
we have plotted the piezometric line, especially the transition point at different times t — 1.6s, t — 1.7 s, 
t = 1.8 s, t = 1.9 s. In figures 25 and 26, the left side to the transition point corresponds to a free surface 
state and the right one to a pressurized except at t = 1.9 s where we can observe two transition points due 
to the pressurized state propagating from the downstream end. The behavior of the two methods are in a 
good agreement and particularly with respect to the localization of transition points. This short and "sharp" 
water hammer test allows us to validate numerically the two approaches for capturing the transition between 
free surface and pressurized flow. 
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Figure 24: Initial state and boundary conditions. 
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Figure 25: Water hammer test case in non uniform closed water pipe. 
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Figure 26: Water hammer test case in non uniform cfosed water pipe. 



5.5 Numerical order of the discretisation error of the kinetic numerical scheme 

In order to obtain a qualitative behavior of the scheme and to compute a "numerical" order of the discreti- 
sation error of the kinetic numerical scheme, we present now a numerical experiment where the steady state 
is not a constant (in x) steady state. The pipe is a circular pipe of diameter 3 m and 100 m long with 
slope 0.001 and Strickler coefficient is Kg — 63.7. The altitude of the upstream end of the pipe is 100 m. 
The upstream total head is kept constant equal to 104 m whereas the downstream water level varies (see 
figure 27). We have compute the "exact" numerical flow for all time from t = s until the stationary state 
is reached at time t = 100 s, by the VFRoe method presented in [2] with a uniform discretisation of 8000 
mesh points. We have then computed the L^ norm of the difference between the piezometric line computed 
by the numerical kinetic scheme for different mesh sizes of the uniform discretisation, Ax, and the "exact" 
numerical solution at time i = 20 s and t — 100 s. 
Other parameters are: 

CFL 

Simulation time (s) 

Sound speed {ms~^) 



0.9, 
100, 
40. 



We present, in figure 28, the piezometric line and the speed of the flow along the pipe at time t — 20 s 

for three different mesh sizes (in fact we prefer to talk of the number of mesh points). The three curves 

representing the piezometric line are very close whereas the coarse mesh does not capture at all the speed 

along the pipe. 

We present, in figure 29, the piezometric line and the speed of the fiow along the pipe at time t = 100 s. 

The three curves representing the piezometric line as well as the speed along the pipe are very close. One 

can see that the stationary speed is not constant along the pipe. 

The numerical order at time i = 20 s, represented in figure 30(a) for different mesh sizes, and the numerical 

order at time t = 100 s, represented in figure 30(b), are almost equal to 1, which was expected since a kinetic 

finite volume scheme is known to be of order 1. 
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Figure 27: Piezometric line at the downstream end of the pipe. 
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Figm^G 28: Piezometric line and speed along the pipe at time i = 20 s. 
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Piezometric line at T = 100 s 
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Figure 29: Piezometric line and speed along the pipe at time t = 100 s. 
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Remark 5.1. Although we do not know if the two numerical schemes that we proposed satisfy the conservative 
in cell entropy (see Equation (6)), every numerical results presented have a very good qualitative behavior. 

6 Conclusion and perspectives 

We have proposed in this work a new manner to extend the numerical kinetic scheme with reflections build 
by Perthame and Simeoni [21], to closed water pipes with varying sections and not only to rectangular closed 
water pipes. This scheme is wet area conservative and under a CFL condition preserves the positivity of the 
wet area. Moreover, by introducing a new interfacial steady profiles which allow the "exact" computations 
of every steady states, we are not restricted in the choice of the x function that is the key point in the kinetic 
scheme and the mathematical and numerical properties of the model (see [21, Theorem 3], Proposition 4.1 
and Theorem 4.1). 
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As a well known feature of general kinetic schemes, we are able to "naturally" deal with flows where a drying 

or flooding zone may be present. This key property was not solved by the previous VFRoe scheme that 

we proposed in [2] without introducing a cut-off function for the wetted area which may causes a loss of 

conservativity. 

The PFS model is numerically solved by a kinetic scheme with reflections using the interfacial upwind of 

all the source terms into the numerical fluxes. Moreover, we have proposed a manner to obtain an exactly 

well-balanced scheme. 

As mentioned in [6, 2] this numerical method reproduces correctly laboratory tests for uniform pipes (Wig- 

gert's test case) and can deal with multiple transition points between the two types of flows. The code to 

code comparison for pressurized flows in uniform pipes has proved the robustness of the method. But due to 

the lack of experimental data for drying and flooding flows, we have only shown the behavior of the piezo- 

metric line which seems reasonable (at less no major difference was observed). For non uniform pipes, the 

two numerical schemes are in a very good agreement even though we are not in possession of experimental 

data. 

We are at the present time interested in the construction of a class of "in cell entropy satisfying" schemes 

consistent with the numerical approximations of hyperbolic systems with source terms. 

The next step is to take into account the air entrainment which may have non negligible effects on the 

behavior of the piezometric head. A first approach has been derived in the case of perfect fluid and perfect 

gas modeled as a bilayer model based on the PFS model [3]. 
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